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ABSTRACT 

A n ew era of directly imaged extrasolar planets has produced a three-planet system (jMarois et alJ 
l2008f ). where the masses of the planets have been estimated by untested cooling models. We point 
out that the nominal circular, face-on orbits of the planets lead to a dynamical instability in ^10^ yr, 
a factor of at least 100 shorter than the estimated age of the star. Reduced planetary masses produce 
stability only for unreasonably small planets (< 2 Mjup). Relaxing the face-on assumption, but still 
requiring circular orbits while fitting the observed positions, makes the instability time even shorter. 
A promising solution is that the inner two planets have a 2:1 commensurability between their periods, 
and they avoid close encounters with each other through this resonance. That the inner resonance 
has lasted until now, in spite of the perturbations of the outer planet, leads to a limit < lOMjup on 
the masses unless the outer two planets are also engaged in a 2:1 mean-motion resonance. In a double 
resonance, which is consistent with the current data, the system could survive until now even if the 
planets have masses of ^ 20 Afjup. Apsidal alignment can further enhance the stability of a mean- 
motion resonant system. A completely different dynamical configuration, with large eccentricities and 
large mutual inclinations among the planets, is possible but finely tuned. 
Subject headings: celestial mechanics — planetary systems — methods: numerical integration 



1. INTRODUCTION 

The method of direct imaging for the discovery 
of extrasolar planets has yielded spectacular first re- 
sults over the last se v eral years dC hauvin e t al.l [200^ 
Lafreniere et al.l 120081: iMarois et all .2008,: .Kalas et al.l 
2008t iLagrange et al.ll2009| ) . Direct imaging is a method 
for discovering planets located far from their host stars, 
an as-yet unexplored region of parameter space, and it 
promises new opportunities to characterize the planets 
using their own radiation. However, because the gravi- 
tational infiuence of directly-imaged planets is not mea- 
sured and the astrometric orbital arcs obtained so far 
are short, determining the planetary masses and orbital 
architectures of these systems is challenging. 

In the newly-discovered planetary system HR 8799 {— 
HD 218396), three planets have been imaged at projected 
separations of 24, 38, and 68 AU from their host star 
(jMarois et al.| [2b08'). The best current estimate of their 
masses is derived from the planetary luminosities, mea- 
sured in the infrared. Because these planets are young 
and massive, they are still radiating prodigiously as they 
contract, cool, and become more gravitationally bound. 
The masses are estimated using untested models of this 
contraction and cooling process. One class of such mod- 
els, the "hot-start" models, provides the largest luminos- 
ity possible at a certain mass and age, given assumptions 
about opacities in the planetary atmosphere. Hot-start 
models have initially extended envelopes and a large en- 
tropy per baryon; even hotter models converge to a com- 
mon track after a few Myr (Baraffc et al. 2002). There- 
fore, for a given age and luminosity, these models should 
provide a lower limit on the mass. For HR 8799, the 
lower-limit masses are 5-11, 7-13, and 7-13 A/jup for plan- 
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ets b, c, and d, respectively, based on a rather uncer- 
tain stellar age of 30-160 Myr^, which is presumably also 
roughly the age of the planets. 

The following simple calculation illustrates why a lower 
mass limit can be inferred from a planet's contraction lu- 
minosity. For HR 8799, the planetary luminosities have 
been measured to be L ~ 10~^Lo, and radii of i? ~ 
1.2i?jiip were derived from the objects' temperatures, 
measured by fitting photometry wi th a variety o f syn- 
thetic spectral energy distributions (jMarois et al.l [2008'). 
Because the objects are cooling, they were more lumi- 
nous in the past, so they have radiated at least i^agc ?J 
4x lO'^^ erg. Their current binding energy, which supplied 
this luminosity, is ~ GM'^R'^ ~ 3 x 10''3(M/Mjup)^ erg, 
where the radius is roughly independent of the mass 
for Jupiter-mass objects. Consequently, M > IMjup. 
Cooling models also take into account that L diminishes 
with time, and thus arrive at a considerably larger mass. 
Whether this larger calculated mass is a robust lower 
limit depend s on th e accuracy of the model. Recently, 
iDupuv et al.l ( j2009l ) measured the dynamical masses for 
a system of brown dwarfs (both of mass w 57Afjup) and 
showed that cooling models overpredict the component 
masses by ^^25%. 

If energy is lost during the process of planet formation, 
then an even larger planet mass would be needed to gen- 
erate the currently observed luminosity. For example, 
in the planetary core-accretion models of iMarlev et al.l 
(j2007| ). considerable luminosity is radiated in the accre- 
tion stream and shock, and that energy is not internalized 

^ This estimate is given bv lMarois eTaLl l l200g) based on four age 
indicators: Galactic space motion, main-sequence fitting, stellar 
pulsations, and the massive debris disk. The first is circumstantial 
but consistent with the quoted ages; the others suggest an age 
< 100 Myr. That the star has reached the main sequence suggests 
that it is at least several 10s of Myr old. 
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by the planet. At the end of formation the planet has less 
gravitational potential energy to later supply its lumi- 
nosity. The integrated luminosity since formation would 
not account for the planet's current binding energy, so 
the mass needed to supply an observed luminosity at a 
given age may be much bigger. 

The HR 8799 system has survived an order of magni- 
tude longer than the primordial gas dis k, which, if typical 
of di s ks of A stars, lasted < 3 Myr ([Hillenbrand et al] 
Il993t [Hernandez et al.l l2005l) . The system has there- 
fore had time to dynamically evolve in the absence of 
gas. Though the planets orbiting HR 8799 are sepa- 
rated by tens of AU, the inferred minimum masses of 
the planets were large enough that their mutual grav- 
itational interactions are important. For example, a 
planet with mass Mp = 10 Afjup orbiting a star of mass 
M, = 1.5 Mq at semi-major axis a = 40 AU domi- 
nates gravitational dynamics within its Hill radius of size 
Rh = a{Mp/3M^y/^ = 5 AU. Because Rh is a large 
fraction of the planetary separation, gravitational inter- 
actions among the planets can substantially modify the 
dynamical evolution of the system. 

In fact, the no minal orbits reported in the discovery pa- 
per (jMarois et al. 2008) are unstable. We integrated the 
Newtonian equations of motion of the proposed system 
using the Bulirs h-Stoer (BS) algorithm of the Mercury 
(|Chamberslll999t ) package (version 6.2), with an accuracy 
parameter of 10~^^. The planets are assigned circular, 
face-on orbits, and we used the nominal masses for all 
four bodies: 7, 10, 10 Mjup for planets b, c, and d, re- 
spectively, and 1.5 Mq for the star^. Figure [T] shows 
the results for the semi-major axis and maximum radial 
excursion of each planet as a function of time. A close 
encounter between planets c and d at 0.298 Myr (i.e., 
they enter within one Hill radius of one another) leads to 
a brief interval of strong scattering which ejects planet b 
at 0.316 Myr (i.e., it reaches > 500 AU with positive en- 
ergy, and is removed from the simulation). Planets c and 
d swap orbits and finish in a stable configuration, with no 
further semi-major axis evolution, but they exhibit a reg- 
ular secular eccentricity cycle with a period of 1.5 Myr. 
This evolution is not unique in its details since the or- 
bital evolution is chaotic. However, qualitatively similar 
evolutions are common for simulated planetary systems 
constructed to match the discovery data: instability usu- 
ally sets in well before the star's age of > 30 Myr. 

The goal of this paper is to determine orbits that are 
consistent with the astrometric data, the inferred plane- 
tary masses, and with dynamical stability over the sys- 
tem's age. Neglecting stability considerations, there is a 
large amount of freedom in fitting orbits to the discovery 
data, because (1) the measured astrometric arcs cover 
only ^2% of the middle orbit and ~1% of the outer or- 
bit, (2) the velocity of the inner planet is almost entirely 
unconstrained, and (3) the line-of-sight positions and ve- 
locities of the planets relative to the star are unknown. A 
priori, two classes of orbital architectures are possible — 
those in which the planets occupy roughly coplanar orbits 
and those with large mutual inclinations. Since plan- 
ets form in disks, it is likely that they initially occupy 
nearly coplanar orbits, and systems that remain stable 

* For the specific initial conditions of this and other integrations 
herein, see Tables m and H 
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Fig. 1. — Semi-major axis, periapse, and apoapse for the three 
planets as a function of time in a numerical integration of the 
nominal model of the system, which has circular, face-on orbits, 
and planetary masses Mj = 7Mj^p, Mc = lOMjup, and = 
lOA/jup. An instability occurs only 3 X 10^ yrs into the integration, 
in which planets b and c suffer a close encounter, suggesting the 
planetary masses or orbits of the nominal model are in error. 

indefinitely are likely to stay roughly coplanar. Alterna- 
tively, the system may not be indefinitely stable. While 
old compared to the lifetime of the protoplanetary disk, 
the current age of the planetary system is probably less 
than one t enth the m ain-sequence lifetime of the star 
(^1.5 Gyr; Ilbenlll96"7l ). Without further analysis, it is 
thus possible that the planets are in the process of scat- 
tering off of one another, currently have large eccentrici- 
ties and mutual inclinations, and will not be stable over 
the lifetime of the star. In fact, current models predict 
that planetary systems undergo periods of strong mutual 
excitation, perhaps generically leadin g to the ejection of 
planets (e.g., [Lev ison et al. 1998, Goldreich et al.ll2^004 
IScharf fc MenolJlIOOgrfVeras et al.ll2009D . 

To explore these possibilities systematically, we take 
the following approach. We start with restrictive as- 
sumptions about the orbital architecture of the system, 
and we then progressively relax those assumptions. At 
each stage, we find parameters that maximize the sta- 
bility, and we finally argue that a resonant configuration 
is most likely for the system to have survived to its cur- 
rent age. In Sp1 16.3( we assume that the orbits of the 
planets are close to coplanar. In fj^l we discuss astro- 
metric constraints on the orbits and show that, although 
the data are consistent with circular and coplanar or- 
bits, the system orientations that fit the data best do not 
generate stable orbits. In ^ we determine what plane- 
tary masses would be needed for circular, coplanar orbits 
to be stable and argue that they are too low given the 
observed luminosities. Having thus ruled out circular, 
coplanar solutions, we next allow the planetary eccen- 
tricities to vary. Since the inner planet's eccentricity is 
unconstrained by the data, we first scan over non-circular 
orbits for the inner planet while keeping the outer two 
planets on circular orbits (SI- The suggestive results of 
this experiment led us to our preferred configuration for 
the planetary system: a mean motion resonance between 
the inner two planets (f}5]). An initial exploration of the 
parameter space of possible resonant orbits shows that 
if the outer two planets are also in a mean motion reso- 
nance, the system could be stable even if the companion 
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TABLE 1 

ASTROMETRIC CONSTRAINTS 



Planet 


[XE, a;jv] (AU) 


s (AU) 


[VE, -wjv] (10-^ AU day-i) 


Vp (10^^ AU day-i) 


b 


[60.16(6), 31.50(6)] 


67.91(6) 


[1.49(0.14), -2.18(0.14)] 


2.64(0.14) 


c 


[-25.90(6), 27.76(6)] 


37.97(6) 


[2.15(0.14), 2.47(0.14)] 


3.27(0.14) 


d 


[-8.45(6), -22.93(6)] 


24.44(6) 


[-3.5(2.7), 0.0(2.7)] 





Note. Sky-projected positions and velocities of e ach planet relative t o the star, found by a 

rectilinear least-squares fit to the astrometry of Table 1 of IMarois et al.l ||2008I1 . Positions are at the 
epoch of 2008 Aug. 12, velocities assume no detectable orbital acceleration. 



masses are twice as large as the nominal masses. In [JH 
we allow all the orbital parameters to vary via a Monte- 
Carlo method. We confirm that mean-motion resonance 
is the most likely reason the planetary system has sur- 
vived. We also find a scattering-type configuration that 
is stable for 30 IMyr, but we argue that it is unlikely. In 
^we discuss our conclusions. 

2. ASTROMETRIC CONSTRAINTS 

In Figure [U we plot the sky-projected position and 
velocity vectors (Table [ij of the three planets, at the 
epoch 2008 Aug. 12, as deter mined by leas t-squares fit 
to the astrometry in Table 1 of llVlarois et al.l ((2008). The 
distance to the star is 39.4 ± 1.1 pc, based on the Hippar- 
chos parallax (^van Lecuwen 2007). We use this nominal 
distance to convert observed angular separations to AU. 
The 3% error thus introduced into the distances and ve- 
locities does not change our qualitative conclusions; in 
^we take this error into account. 

The impression given by Figure [5] is that we are see- 
ing the planetary system face-on, with counter-clockwise, 
nearly circular orbits. This is what we call the "nomi- 
nal model," and we plot the implied orbits and veloc- 
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Fig. 2. — Observed sky-projected positions and velocities for the 
three planets, along with the velocities of face-on, circular orbits 
for M^, = 1.5M0 (the nominal model; model A of Tables [2] and HJl . 
The circles represent 1-cr and 2-cr errors on the measured velocities. 
The inner planet, d, has a barely-detected velocity due to a short 
time baseline. Errors on the positions lie within the circles marking 
the locations of the planets. 



ity vectors for a 1.5 Mq star, also in Figured] If all 
of the orbits are truly face-on and circular, their sky- 
projected separation s = -^/xl + = a, and their sky- 
projected velocity Vp = \/v'e + v'^ = Worb, the orbital 
velocity. Since all of the planets are bound mostly by 
the mass of the star, they should follow circular orbits at 
semi-major axis a with velocities Worb = 27r AU yr~^ 
(Af^/MQ)i/2(a/AU)~i/2. For the outer two planets, 
s and Vp are measured with high precision (Table [1}, 
providing two independent measurements of the stellar 
mass. Given this nominal model, the stellar mass bind- 
ing planet b is Mi^t, = 1.6O±O.17M0 and the stellar mass 
binding planet c is M^^c = 1-38 ± O.12M0. These values 
bracket the value of M^, = 1.47 ± O.3OM0 preferred by 
combini ng parallax, niagnit ude, and spectroscopic infor- 
mation (|Grav fc Kavelfl999l) and are in reasonable agree- 
ment: AM^ EE M^f, - M^^c = 0.22 ± O.2IM0. However, 
there is some tension in the observed velocities. For both 
planets b and c, the observed velocity vector is ~2ct away 
from perpendicular to the separation vector (from the 
star to the planet). The instability reported in the in- 
troduction is, however, the main failing of the nominal 
model. 

To address this failing, we first search for another 
model in which the planets are still coplanar and cir- 
cular, but the system plane is inclined by an angle i to 
the plane of the sky, with an ascending node measured 
East of North, and a to-be-determined consistent mass 
M^. The sky projection changes the magnitudes and di- 
rections of the velocity vectors and the inferred spacings 
of the planets, and taking it into account could lead us to 
infer a wider-spaced, more stable system. We focus only 
on circular and coplanar models in this section, saving 
more complicated direct fits to the data for fj6l The ve- 
locity field on the sky due to this model is: 




and 



n{xE,XN) = {GA'Uy/^ia^ + P^icosiy^y^^^ (3) 

is the mean motion as a function^ of position. 

We solve for the three parameters i, fl, and M^,, as- 
suming that the planets are on non-interacting Keplerian 

^ Here we neglect the few-percent contribution of the planetary 
mass 
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TABLE 2 

Solutions of circular, coplanar systems 



Solution 


(Mq) 


i 


Q 




d.o.f. 


a 


at{AU),Xb 


ac{AU),Xc 


aaiAU),Xa 


A 


= 1.5 


= 0° 


= 0° 


12.63 


6 


0.049 


67.91, 62.36° 


37.97, 316.99° 


24.44, 200.23° 


B 


1.44 


= 0° 


= 0° 


12.19 


5 


0.032 


67.91, 62.36° 


37.97, 316.99° 


24.44, 200.23° 


C 


= 1.5 


21.3° 


151.5° 


9.07 


4 


0.059 


72.89, 62.30° 


38.15, 315.97° 


25.47, 202.23° 


D 


1.86 


33.2° 


145.9° 


5.67 


4 


0.225 


81.04, 61.30° 


38.16, 315.29° 


27.69, 204.92° 


E 


2.28 


41.4° 


143.6° 


2.77 


3 


0.429 


90.02, 60.20° 


38.16, 314.81° 


30.33, 207.30° 



Note. — The symbol denotes values among the parameters AIj, , z, fi that arc held fixed for this solution. 
The inclination i is the angle between the planetary system's orbital angular momentum vector and the vector towards 
the observer, and the ascending node Q is measured East of North (so the position angle of a planet as it passes through 
the plane of the sky, toward the observer, is Q). a is the significance value of being this high using a x^-fest, given a 
certain number of degrees of freedom (d.o.f.), under the null hypothesis that the circular coplanar model with the given 
orientation and stellar mass is correct. For instance, the probability of observing these velocities given the nominal model 
is < 5%. 



orbits, each of which only feels the mass of the central 
star. We calculate values using ve and vn, and their 
associated measurement errors, for all three planets (Ta- 
ble [U 6 data points) — we neglect the errors on xe and 
a; AT, which are too small to affect our results. 

Solutions are reported in Table [2] In model A, which 
is the nominal model, we fix the parameters to their 
nominal values to serve as a baseline. In model B, we 
require face-on orbits, but let M^, float, the result be- 
ing not far from the nominal stellar mass. In model C, 
we fix M^, = 1.5Mq, but let the orientation float. The 
orbits depart from face-on by ~20°, and improves 
a little. In model D, we let all three parameters float, 
but respect the independently measured stellar mass by 
including [{M^/Mq - 1.5)/0.3]2 in x^- In model E, all 
three parameters float with no such mass constraint. The 
orientation-dependence of is shown in Figure [3l and 
the mass-dependence is shown in Figured Interestingly, 
the best fits are for M^, much larger than the nominal 
value 1.5±O.3M0. This is not surprising given the good 
agreement of circular orbits because we are introducing 
line-of-sight offsets and velocities, so a more massive star 
is needed to make such orbits circular. Figure [5] shows 
how the velocity vectors of model D falls into the l-a 
error ellipse for each planet. However, the inner two or- 
bits are closer spaced than the nominal model, and the 
instability is even more rapid: in an integration a close 
encounter occurred between c and d on their second con- 
junction (see Table [3] for initial conditions). 

This integration and all those in §31Hllwere performed 
using the HYBRID integrator of Mercury with a timestep 
of 100 days. Each integration was terminated when any 
two planets passed within 1 Hill radius of each other, 
one was ejected (distance to the star > 500 AU with 
positive energy), or the system lasted 160 Myr. Before 
the onset of close encounters, energy was conserved to 1 
part in ~10^ and angular momentum was conserved to 1 
part in ~10^^. Though we used the HYBRID integrator, 
because the integrations were halted at the first close 
encounter, the integrator's treatment of close encounters 
did not affect our results. In W.ll we verify that after 
a close encounter, at least one planet would be quickly 
ejected. 

Similarly, we fit the best orientation for M^, values 
between I.I-S.OMq, spaced by O.OIMq, and integrated 
those orbits. No three-planet systems generated in this 



contours, for M,= 1.86Mq 
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Fig. 3. — Model ^ function of orbit orientation of cir- 

cular, coplanar models at stellar mass of M* = 1.86M0. Note 
the degeneracy O — > f! -|- 180° , which arises because we are mod- 
eling an unobserved z and i>z (the direction z is away from the 
observer) for each planet, but due to the sky-projection inherent 
in the observations, these values could just as as well be —z and 
— switching the ascending and descending nodes at z = 0. A 
face-on orientation lies at the origin. 
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system orientation (i and Q) of circular, coplanar models. Masses 
above the nominal mass 1.5M0 are preferred. The different crosses, 
each assigned a letter, are different solutions as given in Table [2] 



Multiplanet System HR 8799 



5 




100 50 -50 

E/W separation (AU) 



Fig. 5. — Observed sky-projected positions and velocities for 
tlie three planets, along with the velocity predictions of model D 
(Tables [2] and |4]l of a circular, coplanar system. 

way were stable for more than 1.5 x 10^ yr. Therefore 
we find that more careful fits to the data, under the hy- 
pothesis of circular coplanar orbits, do not simply lead 
to a stable solution. 

3. MUCH LOWER PLANETARY MASSES? 

Before relaxing the assumption that the planets' orbits 
are circular and coplanar, we ask how low the planets' 
masses must be for the nominal orbits to be stable. In- 
tuitively, if their masses are very small, the planets will 
not significantly perturb each other on the timescale of 
30 — 160 Myr. There is a well-developed framework for 
quantifying long-term stability in systems with only two 
planets. In three-body systems, conservation of total 
energy a nd angular mo mentum constrains the possible 
motions (jMarchal fc Bo zis 1982). Applied to a system 
of a star and two planets, we may define Hill stability 
as a constraint that the planet that is initially closer to 
the star stays closer to the star for all time. When the 
criterion for Hill stability is satisfied, a close encounter 
between the planets is prohibited (although escape of 
the outer planet to infinity, or the collision of the inner 
planet with the star, is not forbidden). Qualitatively, 
stability requires that the planets be separated by more 
than a few mutual Hill radii: i?H = ^(oin + aout)e, with 
e = [(Min + Mout)/(3M^)]i/3. Define A as th e planets' 
differe nce in semi- major axes in terms of i?H- iGladmanI 
(flQQl gave the Hill stability criterion as: 

A > A,„. ^ 2V3[l+3V2e-(^ii^l^|±l^) 3-^^6-2+.. 

(4) 

Evaluating these numbers using the nominal system with 
nominal masses 7, 10, 10 Mjup, we have Ad = 2.68 and 
Acrit.cd = 4.03 for the inner two, and Abe — 3.69 and 
Acrit,f)c — 3.98 for the outer two. Apparently both sub- 
systems fail to satisfy the Hill stability criterion. 

We performed numerical simulations to find just how 
small the planets would need to be to remain stable. We 
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Fig. 6. — Time to instability versus scaling of planetary masses. 
In each simulation, the planets are given circular, coplanar orbits 
(model A) , but all of their masses are scaled down by a common 
factor from their nominal values (see Table [4t. Solid vertical lines: 
results for the inner two planets, in the absence of planet b. Dark 
gray region: results fo r all three planets. Light gray region: the 
stellar age as given bv lMarois et al .l (120081). Das hed l ine: t ime at 
which hot-start models of planets from Baraffc e t al.l I I2003I) reach 
luminosities of IO^^'^Lq for masses scaled to -A/nominal = lOMjup. 
Given the system age estimate from [Ma rois et al. (200^, stability 
requires planet masses < 2A^jup. Cooling models allow for stability 
of planets up to ~3.5Afjup at the cost of an uncomfortably low 
system age: ~10 Myr. 



are helped by the long orbital periods and short system 
age (only ^10® dynamical times), which allows suites 
of integrations to be rather inexpensive. First we sur- 
veyed the instability near the nominal orbits ("A"), as 
the search of ^ did not reveal any more stable starting 
points. Let us focus on the inner sub-system (c-d), as it 
is further from stability, and ask the question: by what 
factor must we multiply the nominal masses for stability 
over 30 Myr? In Figure [5] we plot the time to instabil- 
ity — when the first Hill-sphere entry occurs — versus 
this common mass scaling. Vertical lines represent two- 
planet systems consisting of planets c and d on circular, 
face-on orbits. We note that below Mp = 0.33Mnominai, 
where Hill's stability criterion is satisfied, all of the two- 
planet system s last for 1 60 M yr, when the integrations 
were stopped. IGladmanI ()1993D found that if the planets 
initially have small eccentricity (radial excursions com- 
parable or less than a Hill radius) and are not in res- 
onance, then the timescale for instability drops rapidly 
after this boundary (eq. [4]) is crossed. However, insta- 
bility does not appear on timescales relevant for the c-d 
subsystem until Mp ~ O.SMnominai- In a separate suite 
of integrations (not plotted), we found that the nominal 
orbits and masses of the outer pair of planets can be sta- 
ble for 160 Myr. We also plot the instability timescale of 
the three-planet system (dark gray region), with each of 
the three planetary masses scaled by a common factor. 
The masses must be lower than about 1 /5 of the nominal 
masses to remain stable 30 Myr, the lower limit on the 
stellar age (depicted by the light gray stripe). 

When considering three or more planets, there are no 
sharp stability boundaries, but there are well-established 
empirical scaling relations betwee n semi-major axis sep- 
aration and insta. bility timescale ([Chambers et al.lll996l : 
Zhou et all |2007|). A pplying the scaling relation of 
Chatteriee et al.l (|2008l appendix A, fit 1) to the HR 8799 



6 



Fabrycky and Murray-Clay 



Sd pc 















0.6 0.8 1.0 1.2 1.4 1.6 

0^(24.44 AU) 

Fig. 7. — Instability time of coplanar systems, with planets b 
and c on circular face-on orbits and planet d on a non-circular orbit, 
with masses of 5, 7, and 7 Mj^p for b, c, and d respectively (see 
Table|4]l. Lines: based on planets c and d (in the absence of planet 
b) from the nominal model, but choosing a non-zero eccentricity 
for planet d to satisfy its currently-observed distance from the star. 
Dark gray region: same as before, but in the presence of planet b 
with its nominal parameters . Light gray region: the stellar age as 
given bv lMarois et al.l II2008I ). A few three-planet systems last the 
stellar age, and these correspond to a 2:1 mean motion resonance 
between planets c and d. Current separations between the star and 
the planets are labeled s^j and Sc- 

system implies A > 4.4 if the system is to remain stable 
> 30 Myr. Let us assume the instability between the in- 
ner two planets is dominant, so this limit applies to Ad, 
then the masses must be < 1/4 of the nominal masses, 
in good agreement with the non-resonant systems of Fig- 
ure[Sl We note, however, that none of the published scal- 
ing relations extend to planetary-to-stellar mass ratios of 
5 X 10~^, nor do they strictly apply if adjacent planets 
are unequally spaced in Hill radii {A^d 7^ A{,c), both of 
which are relevant for HR 8799. For masses > 1/4 of the 
nominal mas ses, our results give a lo nger lifetime than 
the scaling of IChatteriee et all ()2008[ ). sometimes orders 
of magnitude longer. Regardless, we confirm that insta- 
bility can occur even if the sub-systems are initially Hill 
stable. In circular, face-on orbits, the implied upper lim- 
its of masses — 1.5, 2, and 2 Mjup — are incompatible 
with any cooling model at ages greater than 30 Myr, even 
extreme hot-start models. 

4. NON-CIRCULAR INNER ORBIT? 

In the previous section, we found that face-on, circu- 
lar orbits, consistent with the astrometric constraints, 
could only be stable if the planetary masses were im- 
plausibly low. In this section, we choose the lowest plan- 
etary masses that are compatible with hot-start models, 
and we choose a non-circular orbit for planet d (its or- 
bit is currently unconstrained by observations). We seek 
systems that remain stable until the lower limit on the 
stellar age of 30 Myr. 

We first simulate the inner two planets, each of 7Mjup, 
in the absence of planet b. They are given copla- 
nar orbits, with planet c on a circular orbit at Oc = 
Sc — 37.97 AU. The initial longitudinal separation is 
given by the observed positions, assuming face-on orbits 
(Ac — Ad « 117°). We scan over a grid of semi-major 
axes for the inner planet. For ad < Sd — 24.44 AU, 
Cd is chosen so that apastron is at 24.44 AU, and for 
Qd > 24.44 AU, Cd is chosen so that periastron is at 
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numerical integration with = 0.95 X 24.44AU = 23.22AU from 
Fig. [7] (see also Table [4]|. The rotating coordinates are centered 
on the star with planet c on the positive horizontal axis. The 
circle is a distance from the star of 24.44 AU. When the inner 
planet lags the middle planet by ~ 117° (its current position), its 
distance from the star is ~ 24.44 AU, but when it reaches the 
middle planet's longitude, it is always closer to the star. Planets 
are labeled near their currently-observed positions relative to one 
another. The 2;1 mean motion resonance protects the two planets 
from close encounters and adds coherence to the long-term transfer 
of energy and angular momentum during encounters. 

24.44 AU (see Table |4] for how the initial conditions are 
generated). These choices maximize the chance that the 
two-planet system will be stable, while matching the con- 
straint of the currently-observed separations from the 
star. We plot the instability times in Figure [7] as vertical 
lines. We repeat this calculation with planet b present 
with its nominal orbital elements (see Table H]) and with 
mass 5 Mjup, and plot those instability timescales in Fig- 
ure [7] as a gray region. 

We observe that a very narrow range of is compat- 
ible with both the observed astrometry of the planets 
and with dynamical stability. The presence of the third 
planet narrows this range still further. The center of this 
range corresponds with the 2:1 mean motion resonance 
between planets c and d. (The position is offset from the 
location = (l/2)^/'^ac because the large mass ratios 
induce fast precession.) In Figure [H] we show how this 
resonance protects the planets from close encounters. 

We ran identical simulations with planetary masses 
of 7, 10, and 10 Mjup, and found qualitatively similar 
results, except the most stable three-planet simulation 
lasted only 10 Myr. In the next section we examine 
this resonant protection mechanism and find initial con- 
ditions that produce acceptably long survival times even 
for these and even higher masses. 

5. MEAN MOTION RESONANCE 

Inspired by the fact that Figure [7| shows a region of 
greater stability in the vicinity of the 2:1 resonance be- 
tween c and d, we search for a face-on system near 
the center of the resonance. We use the BS integrator 
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throughout this section. We find a sohition in the ab- 
sence of planet b in which the resonance angle, 



(5) 



librates with small amplitude around 0°. When evalu- 
ating resonant angles, we compute the orbital elements 
with astrocentric coordinates. Resonance requires that 
ad is low enough for planet d's period to be commensu- 
rate with planet c's, and the observations require that 
Cd is high enough for planet d to reach its currently ob- 
served separation from the star. Currently, we observe 
\c-\d~ 117°, so 0d sa 0° implies Xd-voa ~ 126°. Thus 
we find small libration is compatible with planet d be- 
ing closer to apoastron than to periastron at the current 
time, in which case the velocity should be smaller than 
that of a circular orbit at the same distance. Integration 
with these initial conditions for planets c and d, in the 
absence of b, indeed shows libration and long-term sta- 
bility (at least 160 Myr), for initial ad and values that 
place the planets in resonance (e.g., the system labeled 
"two-planet resonance" in Table [3]). 

In such solutions, the resonance angle for planet c usu- 
ally does not librate. The resonance involves only the 
eccentricity of planet d. When planet b is added, plan- 
ets b and c excite each other's eccentricities and cause 
the libration amplitude of 4>d to fluctuate. Sometimes 
these excited eccentricities cause an encounter between 
b and c; sometimes the loss of libration in (pd allows an 
encounter between c and d. In Figure [9] we show an ex- 
ample of this instability, where planets b and c start in 
their nominal orbits, ad = 23.32 AU, — 0.09, 0^ = 0°, 
and all bodies have their nominal masses (see Table [3]). 
Panel (a) shows the range of motion of each orbit versus 
time, panel (b) shows the resonance angle versus time, 
and the bottom panels show brief segments (3 x 10^ yr, 
at times labeled above panel b) of the motion of the res- 
onance angle through phase space. Over such brief inter- 
vals, the libration amplitude holds rather steady, except 
at the very end of the integration. In this example, the 
instability causes an encounter between planets c and d 
at 35.6 Myr. 

Compared to the non-resonant cases, this system 
showed considerable longevity: it lasts long enough to 
be a plausible model for the observed system. We have 
found a way to calm the strongest interactions, those that 
cause instability after a few thousand orbits: a resonance 
between planets c and d that protects them from close 
encounters. This resonance protects the system until the 
somewhat longer timescale interactions between b and 
c cause an instability. But those interactions can also 
be suppressed by postulating yet another resonance. We 
integrated the nominal masses with initial conditions as 
above except ad = 23.42 AU instead of 23.32 AU (Ta- 
ble [3]) . The resulting system showed resonance protec- 
tion between planets b and c. The 2:1 resonance is ac- 
tive, which this is possible far from its nominal location 
because the pericenters are precessing on nearly orbital 
timescales. In Figure fTOl we show this system lasting for 
160 Myr. In this example, the resonance angle (pd is li- 
brating with small amplitude the whole time (panels b 
and c), and the resonance angle 0c, out — 2Af, — Ac — njc 
spends more time near 0° (panels d and e), indicating 
the system is protected by both resonances. Even af- 
ter 160 Myr of evolution, we have verified that there are 
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Fig. 9. — A simulation of the nominal masses, which is initially 
protected from close encounters by the 2:1 resonance between plan- 
ets c and d, but it is destroyed after 35.6 Myr due to interactions 
between planets b and c (see Table[3]for initial conditions). (Panel 
a) Semi-major axis, periapse, and apoapse for the three planets as 
a function of time; (panel b) dots: resonance angle every ~ 10^ yr 
(libration is rapid, on nearly orbital timescales, and is not well- 
sampled) and lines: its running envelope, as a function of time; 
(panels tl-t4) phase plot of the resonance angle, over short dura- 
tions, as labeled in panel b. 



epochs at which this solution fits the astrometric data 
of Table [TJ We found print-outs for which a rotation in 
the plane of the sky matched the simulated to the ob- 
served positions within a fractional error of 1% (more 
print-outs would likely find a closer match), and then we 
calculated based on the velocities of Table [T] The 
resulting — 11-4 was both acceptable and quite com- 
petitive with the models of Sj2l 

The next step is to understand how these resonances 
protect the system as a function of planetary mass. 
For instance. Figure [7] shows four integrations in which 
the resonance allows planets of masses M{, — 5, and 
Mc = Md = 7 Mjup to be stable for 30 Myr, which is 
consistent with the observed system. But can the system 
survive at the nominal masses with only one resonance? 
How high can the masses go, in the double resonance? 
In Figure [11] we plot the time to instability for a wide 
range of planetary mass scalings. We use initial condi- 
tions corresponding to the nominal face-on, circular or- 
bits (non-resonant), the initial conditions for Figure [9] 
(singly resonant), and parameters chosen to maximize 
stability of the double resonance for massive planets. All 
are listed in Table [31 Because the resonant locations shift 
with increasing planetary mass, the ideal orbital parame- 
ters for stable resonance depend on the masses. In a suite 
of integrations we slightly vary the initial conditions (see 
Table [3]) to sample the chaotic outcomes. 
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Fig. 10. — A simulation of the nominal masses, lasting for 
160 Myr with no signs of imminent instability, due to a double 
resonance (see Table [3] for initial conditions). Panels are as in 
Figure |9] The resonance angles are defined as 0^ = 2 Ac — A^j — ro^j 
and (pcout = 2A5 — Ac — roc. 



We find that systems with the nominal masses rarely 
survive 30 Myr with a single resonance, but can easily 
survive at least 160 Myr with a double resonance. In 
fact, our integrations show that a doubly-resonant sys- 
tem can be stable for 160 Myr, even for planetary masses 
a factor of two larger than the nominal values. That is, 
if this doubly- resonant configuration is correct, the plan- 
ets could even have the masses of brown dwarfs. We 
find it remarkable that a double 2:1 resonance can allow 
planetary masses an order of magnitude larger than the 
'^2Mjup allowed in a stable, non-resonant system. 

In some of these integrations, we have found the three- 
body Laplace resonance, with angle (f)L = Ad — 3Ac -t- 2Ab, 
librating temporarily (see also ^16. 21 below). Such solu- 
tions are also consistent with the astrometric data. The 
Laplace angle w a.s first observed to librate i n the satellites 
of Jupiter (e.g., iMurrav fc DermottI [19991) . Besides HR 
8799, two other extrasolar planetary systems have been 
proposed to inhabit the Laplace resonance. Extra peaks 
in the perio dogram of radial ve locity residuals of the GJ 
876 system (jRivera et al.ll2005[ ) and the HD 82943 sys- 
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Fig. 11. — Time to instability of three-planet systems as a 
function of a common mass scaling for the planets relative to their 
nominal masses. See Table \3\ for initial conditions. To test the 
sensitivity of our results to small changes in initial conditions, for 
each mass scaling, we calculate the time to instability for the stated 
orbital parameters and for five additional sets of orbital parame- 
ters generated as follows. For each orbital element, we draw a 
random number from a normal distribution with mean and a 
standard deviation of 1. We then multiply the result by a scaling 
factor and add it to the initial value of that element. The scaling 
factors are 10~* AU for the semi-major axis, 10~* for the eccen- 
tricity, and 0.01° for the inclination, ascending node, longitude of 
pericenter, and mean anomaly. Plusses: Face-on, circular orbits; 
Diamonds: Orbits in which the 2:1 resonance between planets c 
and d is active initially; Crosses: Orbits in which both the 2:1 
resonance between planets c and d and the 2:1 resonance between 
planets b and c are ac t ive in itially; Gray region: the stellar age as 
given bv lMarois eTaP 1(200^ ). 



tem (jBeauge et al.l 120081) could correspond to planets in 
the Laplace resonance with the known planets. Although 
each of these three extrasolar systems taken separately is 
merely suggestive of 4:2:1 and Laplace resonances, taken 
together they are quite intriguing. They may point to 
a new area of research in multiplanet systems that has 
been explored rather little so far, both theoretically and 
observationally. 

6. MONTE-CARLO SEARCH 

In the integrations so far, we have systematically var- 
ied a few parameters, concluding that a mean motion 
resonance is a promising solution to stabilize the system. 
Now we seek alternatives by allowing all the other or- 
bital parameters to vary. The objective is to survey what 
orbits are allowed when the age of the system and the 
planetary masses are presumed to be robust. To be con- 
servative, we adopt the youngest system age (30 Myr), 
corresponding to the least massive planets (5, 7, 7) Mjup, 
as in SjH We select all the other variables with a Monte 
Carlo method. We draw the stellar mass from a nor- 
mal distribution with mean lAlMp) and standard de- 
viation 0.30Mo ()Grav fc Kavelll999f ). and we draw the 
system distance from a normal distributi on with mean 
39.4 pc and standard deviation 1.1 pc (jvan LeeuwerJ 
|2007| ). The planetary sky-projected positions and veloc- 
ities are drawn from normal distributions according to 
the observed parameters of Table [T] Note that these pa- 
rameters are derived from the discovery observations of 
iMarois et all (|2008) only; in H6.6l we check which systems 
are consistent with the important precovery observation 
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of planet b bv lMetchev et alj (|2009[ ). What remains is to 
draw z and Vz for each planet, which are not constrained 
by the observations; we make that choice in various ways 
in the following subsections. 

In this section, we follow some planets with very high 
eccentricities and in some cases integrate through close 
approaches between planets. We use the BS integrator 
as before, and we follow the integration until one planet 
is ejected or collides with the star. Over 30 Myr, energy 
is typically conserved to better than one part in 10^, and 
angular momentum is typically conserved to better than 
one part in 10^. 

6.1. From crossing orbits to ejection 

First, we wish to verify that once planets' orbits begin 
crossing, at least one of them is ejected in a timescale 
much shorter than the age of the system. To do so, 
note that the expression for orbital energy of a single 
planet around a star is a monotonically increasing func- 
tion of \z\ or \vz\- Therefore, selecting non-zero values for 
those parameters will lead to a planet that is less bound 
than in the case of a face-on orbit, li — 0, there 
is a maximum value of \z\, called jzlmax, that permits 
a bound orbit. If z is given, then there is a maximum 
value of \vz\, called |wz|max, that permits a bound orbit. 
We wish to find how long it takes for planets on cross- 
ing orbits that are not marginally bound to be ejected 
by each other, so we first selected z from a distribution 
uniform in the interval [— |2;|max/3, l^lmax/S], after which 
we selected Vz from a distribution uniform in the inter- 
val [— l^zlmax/S, l^zlmax/S]. This choicc of distribution 
has the advantage of being connected to the observables 
(actually, complementary to them) , being easy to imple- 
ment, and being tuned to answer the question. It has the 
disadvantage of not corresponding to a simple distribu- 
tion in orbital element space; nevertheless, a very wide 
range of orbital elements are sampled. 

We ran 1530 systems generated in this way, integrating 
to an ejection of one component (defined as reaching an 
orbital distance > 500 AU with positive energy). The 
median time to ejection was 0.22 Myr, and the max- 
imum time was 7.7 Myr. These timescales are longer 
than the ~0.02 Myr scattering phase of Figure [1] due to 
significant mutual inclinations. In any case, the scatter- 
ing phase will not contribute significant longevity to the 
system, and we are justified in stopping integrations at 
the first clo s e app roach in other sections of this paper. 
IVeras et"al1 (|2009t ) have also reached this conclusion for 
the planets of HR 8799. 

6.2. Moderately eccentric, coplanar planets 

Next, we extend the analysis of fJH to the case in 
which all three planets have non-zero eccentricities. For 
simplicity, and acknowledging that the planets probably 
formed in a common, fiattened disk, we first investigate 
coplanar systems. There are several steps to generating 
a (z, Vz) pair for each planet: 

1. draw stellar mass, distance, and planetary sky- 
projected positions and velocities, as described 
above; 

2. draw a vector uniformly from the unit sphere, 
which serves as the direction of all the planets' an- 
gular momenta; 



3. compute z and Vz for each planet, consistent with 
the already-chosen spatial variables; 

4. discard the system if x, a number drawn from a uni- 
form distribution in [0, 1], is greater than C/Cml, 
where 

II e,cxp[-{e,/<jf/2] (6) 

j — b,c,d 

and 

/:ML = (aexp[-l/2])3 (7) 

with a = 0.05. 

If a system is discarded at steps 3-4, the process be- 
gins anew with step 1. Step 4 is a technique called re- 
jection sampling, and its purpose is to impose a prior 
distribution on the selected orbital elements et, ec, and 
Cd- We sought planetary orbits with low to moderate 
eccentricity, using the Rayleigh distribution (eq. [B]), as 
recommended by recent work on the generation of ec- 
centricities b y plan etary perturbations (Zhou et al. 200^ 
lJuric fc Tremaineii2008l ) . We chose a consistent with the 
dynamically "inactive" population of lJuric fc Tremaind 
(2008). 

We ran 16, 581 systems generated this way, until a close 
approach or 30 Myr elapsed. The median time until close 
approach was 3,100 yr, and only 49 survived 30 Myr. 
Of the surviving systems, we checked for 2:1 resonances, 
two resonant arguments for the inner pair and two res- 
onant arguments for the outer pair. We considered the 
resonance to be dynamically significant ii h = e cos (f> 
had a non-zero average value (as in Figure ITOl panel e): 
\{h)\ > 2.5y/ {It? /n), where the averages were performed 
over n (^100) printouts of astrocentric orbital elements. 
This criterion is considerably looser than traditional def- 
initions of being "in a resonance," either libration of a 
resonant angle or lying interior to a separatrix in phase 
space. Nevertheless, this criterion indicates (i) protec- 
tion against close encounters in the sense of Figure [5] and 
(ii) enhanced coherency to energy and angular momen- 
tum transfers during conjunction, as conjunctions occur 
at preferential phases of the orbit. All 49 survivors had 
at least one of the four angles fulfilling this criterion: for 
26 only the inner pair were engaged in the resonance, for 
2 only the outer pair were engaged in the resonance, for 
21 both resonances were active. We expect that most 
of these 49 survivors will eventually be disrupted by the 
perturbations of planet b; in no case was the motion as 
periodic as in Figure[TUl Even systems stable for 160 Myr 
may become unstable over the main-sequence lifetime of 
the star (Gozdziewski & Migaszewski 2009). Although 
21 integrations displayed the 4:2:1 double resonance, in 
only one case did the Laplace angle 4>l librate the entire 
time (about 180°), and that system shows the strongest 
inner and outer resonances of the entire set. It is un- 
clear if or how libration of 4>l enhances stability, over 
and above each pair of 2:1 resonances. 

One surviving system showed only a very weak in- 
ner 2:1 mean motion resonance. We examined the first 
1.6 Myr of this integration in detail, finding that Pc/Pd 
fiuctuated in the range 2.30 — 2.45 and Pb/Pc fluctuated 
in the range 2.7 — 3.0. With thousands of print-outs 
we determined that h associated with (j)d had a non-zero 
average of 4-0.012. However, its large range —0.078 to 
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Fig. 12. — A simulation lasting for 109 Myr with weak mean 
motion resonance, but with strong apsidal locking (see Table |4] for 
initial conditions). Small crosses correspond to the apoapse posi- 
tion of each planet, at snapshots spaced by 1.1 Myr, in a frame 
centered on the star with planet c's apoapse on the negative hori- 
zontal axis. The straight lines from the star are the current lines of 
apoapse, the ellipses are the current orbits of the planets, and the 
black dots with white rims are the current positions of the planets. 
All three apsidal lines advance with the same average period of 
1.5 Myr. The azimuthal distribution of the crosses indicates the 
variation in Aro, the radial distribution of the crosses indicates 
the variation in eccentricity (until the very end, the semi-major 
axis change is negligible). For each pair taken separately, the inner 
planet's apocenter tends to be co-located with the outer planet's 
apocenter, so close approaches are forbidden. Note that at the cur- 
rent time planets b and c are near pericenter and planet d is near 
apocenter, enhancing their dynamically-packed appearance. 

-1-0.105 suggests weaker protection by the 2:1 resonance 
than that enjoyed by the other stable systems. The outer 
two planets occupied the 3:1 mean-motion resonance as- 
sociated with the angle i = 3Ah — Ac — 2n7c at a simi- 
larly weak level. Most strikingly, all three planets of this 
system maintained apsidal alignment with one another 
throughout the integration: their relative apsidal angle 
Ataj librated around 0°. Apsidal locking apparently pro- 
vides additional protection against close approaches: an 
inner planet only comes to apocenter at the same spatial 
location where an outer planet is at apocenter, so the 
two bodies do not come too close together. This system^ 
is illustrated in Figure I12[ which gives a pictorial repre- 
sentation of the apsidal protection mechanism. We also 
searched all of the stable systems for a tendency towards 
apsidal alignment, as quantified by |(cos Atu)! being non- 
zero in the same way as h above. We found that apsi- 
dal alignment was also common for systems with strong 
mean-motion resonance; in the set of 49 survivors it oc- 
curred 3 times between the inner planets only, 12 times 
between the outer planets only, and 10 times among all 
three planets. 

6.3. Arbitrarily eccentric, coplanar planets 

We repeated the procedure of §6.21 without step 4, im- 
posing no prior on the eccentricities. However, to ensure 
that close encounters were not already happening, we 

® Initial conditions are given in Table l4l 



only accepted each system if 0^(1 + Cd) < 0.85ac(l — Cc) 
and ac(l -I- Cc) < 0.85ab(l - Cf,). We integrated 25,280 
systems, of which only 5 lasted 30 Myr, and the median 
time to a close approach was 37,000 yr. Of the 5 sur- 
vivors, one had the 2:1 mean-motion resonance active 
and secular alignment between the inner two planets. In 
the other 4, all three planets tended towards alignment, 
and two of these had the 2:1 mean- motion resonance ac- 
tive between the outer two planets. In 3 cases, the inner 
planet had Cd > 0.95 and a current position near apoc- 
enter, and the apsidal alignment between d and c was 
rather tight: \wc — tn^j < 45°. These systems may cor- 
respond to the non-linear secula r resonance identified by 
iMichtchenko fc MalhoTral ()2004[ ). We consider it unlikely 
that configurations with strong apsidal alignment but no 
mean-motion resonance correspond to the true system, 
as explained below (! j6.6p . 

If all three planets orbit in the same plane, one may 
wonder whether the debris disk and the stellar equa- 
tor share it as well. From the observed positions and 
velocities, we can use this Monte Carlo study to con- 
strain the orientation of the planets. Coplanar systems 
that fit the observed positions and velocities with arbi- 
trary eccentricities, but non-crossing orbits, have line- 
of-sight inclinations less than 45°. The subset of stable 
systems obey this same limit. Note that this is not sub- 
stantially different from the limit for circular orbits, cf. 
Figure [3] Constraints on the stellar spin orientation — an 
expected rotational veloci ty v and a measured v sin i — led 
iReidemeister et all ()2009 ^ to derive a stellar inclination 
of 13° — 30°, consistent with this limit. 

6.4. Circular, non-coplanar orbits 

Next, we investigate whether moderately non-coplanar 
orbits are substantially more stable, even in the absence 
of any resonance. If not, we expect the conclusion that a 
resonance is needed applies to any roughly coplanar sys- 
tems, not just strictly coplanar ones. The following series 
of integrations is not intended to match the currently ob- 
served positions, but to be a parametric study of mutual 
inclination. Relative to the nominal case (A), but now 
with lower-limit masses, we varied the initial inclination 
and node of planets b and d (see initial conditions in 
Table |3|). Both planets b and d are given the same in- 
clination (relative to c, which always starts at ic — 0). 
However, they are given different nodes, to sample the 
same inclination 36 times, so that the spread of chaotic 
outcomes are represented. The systems were integrated 
until a close approach or 160 Myr. 

In Figure [13] we plot the resulting times of instabil- 
ity. At low inclination the spread of these times is sev- 
eral orders of magnitude. Varying the initial orientation 
angles, which also control the initial longitude of each 
planet, causes some systems to have a close approach 
within a few tens of orbits and causes other systems to 
last > 1 Myr due to the protection afforded by the 2:1 
resonance between planets d and c. As the initial inclina- 
tion increases, one might expect Hill-sphere encounters 
will be delayed, as the motion out of the plane exceeds 
a Hill radius at i ~ 6° for these masses. Nevertheless, 
we observe that the median instability time modestly de- 
creases until 40°, is constant between 50° and 140°, and 
increases dramatically from 150° to 180°, perfectly retro- 
grade. The shorter instability timescale for substantially 
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Fig. 13. — Inclination dependence of instability times. Planet 
b and planet d are both initially inclined relative to planet c by 
the indicated inclination, and a 6 X 6 grid of initial conditions for 
Ut and in [0°, 60°, ... , 300°] was performed— see Table [3] 
Increased mutual inclinations do not increase stability times for 
prograde orbits. 

non-planar systems is likely due to the iKozail ()1962l) ef- 
fect, whicli causes inclination to decrease and eccentricity 
to increase on a secular timescale of ~10^ years. 

We conclude that we can likely extend our conclu- 
sions from strictly coplanar systems to systems with 
large prograde inclinations. Systems with adjacent plan- 
ets orbiting in the opposite sense can avoid instability; 
i.e., retrograde systems are inherently more stable than 
prograde systems, as has been shown in ot her contexts 
(|Nesvornv et"aII[2rM IcTavon fc BoisI [20081) . No mean- 
motion resonance or apsidal locking appears to be active 
in protecting retrograde systems from instability, but the 
brief timescales of conjunction may be responsible. We 
do not expect a retrograde configuration for planet d, nor 
is it consistent with the data (see 

6.5. Arbitrary orbits 

In this subsection, we remove all restrictions on in- 
dividual eccentricity and orbital orientation, to have a 
completely data-driven sampling of orbits. We expect 
that not many systems generated this way are realistic, 
as planet d's orbit is so poorly constrained. Neverthe- 
less, after drawing the stellar and observable planetary 
parameters, we drew each z uniformly from the interval 
[— 1 2^1 max, 1 2^ I max] and then drew each uniformly from 
the interval [— Iwzlmax, I'^zlmax]- We rejected the resulting 
system if either: 

1. any of the planets initially had positive energy (de- 
spite trying to avoid this case by construction of 
the osculating orbital elements); or 

2. planetary orbits crossed, i.e., adi^+ed) > ac{l — ec) 
or ac(l -I- Be) > ab{l - eb). 

We integrated 3010 systems until the first ejection (not 
stopping at close approaches), or until 30 Myr elapsed. 
Of these, 13 systems survived the entire time. Three 
cases had inner apsidal alignment, and another two cases 
had outer apsidal alignment. A common characteristic 
is that planet d has a moderate-to-large eccentricity and 
comes to apocenter at its currently-observed position. 
Thus the period ratio Pc/ Pd can be quite large, and the 
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Fig. 14. — A non-coplanar system lasting 110 Myr before a 
close approach, shown at the current epoch. See Table|4]for initial 
conditions. The front view (left panel, observer's prospective) and 
right-side view (right panel) are given. The parts of the orbit on 
the observer's side of the sky plane are solid, parts on the far side of 
the sky plane are dashed. The apocenter lines are drawn, showing 
that planet d is currently near apocenter, and planets c and d are 
currently near pericenter, which is usually required of models that 
do not rely on protection from a mean-motion resonance. 

system rather hierarchical. In fact, we found that using 
the f ull formulation of Hill stability ([Marchal fc BozisI 
|1982|) , the Star-d-c sub-system (neglecting planet b) was 
usually Hill stable if planet d is prograde with respect to 
the others. With only small perturbations from planet 
b, the system remains stable for 30 Myr, with the inner 
subsystem fulfilling Hill's stability criterion for most of 
that time. In contrast, as mentioned above, no partic- 
ular protection mechanism was apparent for the retro- 
grade systems. The greater stability of retrograde or- 
bits is not seen in Hill's stability calculations, just as the 
greater stability of retrograde satellites is not reflected 
in the usual Jacobi constant in the circular restricted 
three-body problem (jHamilton fc Krivovlll997[ ). 

A common attribute is that in all 13 stable systems, 
planet b is in a very wide orbit, but it comes to pericenter 
at its currently-observed location. One way to quantify 
this is to note that, in all cases, the mean anomaly of 
b at the present epoch is within 10° of 0°, which would 
happen only 1/18 of the time for randomly phased or- 
bits. In such an orbit, it perturbs the inner two planets 
minimally, contributing to the system stability. The or- 
bits of one of these 13 systems at the current epoch is 
displayed in Figure [HI 

We regard these systems, with very large mutual incli- 
nations and eccentricities, to be a possible, but not plau- 
sible, explanation of the observational data. The outer 
orbit being very near periastron seems finely tuned, be- 
cause it does not spend much time there. It is as if it 
is swooping in from sever al hundred AU, just in time to 
have its picture taken bv iMarois et all ()2008[ ). We also 
note that orbits like those of Figure [T4| are currently con- 
sistent with circular, coplanar orbits, but they do not ful- 
fill this condition at most other orbital phases. Since the 
observed velocities are consistent with circular and copla- 
nar, we can quantify the likelihood that the true system 
actually has very non-coplanar orbits, as follows. We 
produced outputs for the system displayed in Figure [T4| 
at every ^ 3500 years of a 10^ year integration, dur- 
ing which the semi-major axes and eccentricities showed 
no qualitative long-term changes. For each of the out- 
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puts, we found the of the best-fitting circular, copla- 
nar model, minimizing over (M*,i,r2) as in fj2l with 
in the range [1,2] x Mq with no penalty for being far 
from the nominal value. We assign the observed error 
bars on velocities to the sky-projected velocities in the 
simulation, then compare them to the best-fitting ve- 
locities that come from a circular, coplanar hypothesis. 
With the current snapshot (Figure [14]), the system fits 
a circular, coplanar model with = 2.18. In only 639 
of the 28,996 snapshots was better than this value. 
So we conclude that the current phase of this particu- 
lar non-coplanar system is fine-tuned to the ^-^2.2% level. 
We thus prefer models that actually are close to circular 
and coplanar. 

Another consideration disfavors these solutions with 
large mutual inclination and large eccentricities. Fits 
to the infrared spectral energy distribution of HR 8799 
suggest a population of colliding, dust-forming bodies 
with a semimajor axis of ^100 AU, though thi s distance 
is still uncertain dW iUiam s fc Andrews! l2006t ISu et al.l 
l2009t iReidemeister et al. 2009| ). It is unlikelv that planet 
b's orbit actually crosses this belt of bodies, which may 
constrain e\, to less than a few tenths. This constraint 
will be observationally accessible in the near future. 

6.6. Comparison to New Data 

The foregoing analysis was based solely on the astro- 
metric measurements reported in the discovery paper. 
Between the original submission of this work and now, 
several new measurements were reported based on careful 
analysis of previously-collected images (iLafreniere et al.l 



l2009tlFukagawa et ai.ll2009l : lMetchev et al.l2009D. In par 

ticular, the measurement bv lMetchev et al.l ()2009D of po- 
sitions of planet d over a one- year baseline determines its 
sky-velocity to be 

[ve,vn] = [-4.5 ±0.8, 0.7 ±0.8] x 10"^AUday~\ (8) 
These values come from a re gression of po s ition versus 
time, in combination with the lMarois et al.l ()2008[ ) data, 
as in Table [T] For the purposes of this section, that 
value eliminates some of the previously possible dynam- 
ical configurations, as follows. In Figure [T51 we plot the 
l-cr and 2-cr contours of the measured velocity of planet 
b, and overlay the orbits from the preceding sections. 
We find that solutions with a retrograde planet d, and 
those for which planet d is at apastron of an very eccen- 
tric orbit, are now ruled out. However, the 2:1 resonant 
solutions are still quite consistent with the new data, 
including the weakly- resonant, apsidally-locked system 
(Fig 1121) and some of the very non-coplanar (yet non- 
resonant) solutions (e.g.. Fig [T4l) . Although these latter 
solutions fit the data, recall that we are seeing these sys- 
tems at a special time, so we doubt they correspond to 
the true system. 

7. DISCUSSION 

We have investigated the orbital stability of the newly- 
imaged planetary system HR 8799. The nominal orbital 
model and masses are not stable. In fact, no model with 
circular, coplanar orbits that also fits the astrometry well 
is stable, regardless of the inclination and orientation of 
the system on the sky. 

To overcome this problem by reducing the planetary 
masses, values < 2 Mjup are required. This can hap- 
pen if the cooling models under-predict the luminosity. 




6 4 2 -2 -4 -6 
Ve [10"^ AU day"'] 

Fig. 15. — Sky-projected velocity of planet d from vario us sta- 
b le systems th at fit the discovery data l IMarois et al.ll2008l) . from 
§3621 16.31 16.51 The shaded regions are the l-cr (dark gray) and 
2-a (light gray ) regions allowed by the precovery observation by 
IMetchev etaU |(2009). The open circle is the nominal (circular, 
face-on, = I.5M0) solution. Small dots are the solutions with 
2:1 mean-motion resonance between either adjacent pair or a 4:2;1 
double resonance. Crosses are coplanar solutions for which no 2:1 
mean-motion resonance was active, but apsidal alignment was pre- 
served among the inner two planets. Aster isks are solutions with 
arbitrary eccentricity and inclination from i|6.5l 

though that is difficult to understand, as even hot-start 
models cannot produce the observed luminosity at such 
low masses (see SJT]). Such masses would be plausible if 
the system is considerably younger than expected, yet 
the s tar has reached the main sequence (jMarois et al.l 
[2OOI . 

Our favored solution is that a 2:1 resonance between 
the inner two planets preserves stability. Assuming that 
the inner pair of planets are in resonance, two qualita- 
tively different configurations are possible: 

• The outer two planets are not in resonance. This 
configuration remains stable in the perturbing pres- 
ence of planet b only if the planetary masses are 
< 10 Mjup (Fig. [11]). (It is also possible, given a 
less likely system orientation, that only the outer 
resonance is active. This configuration leads to a 
similar mass limit.) 

• The outer pair of planets are also in 2:1 resonance. 
This solution fits all the current data for the sys- 
tem. At the nominal masses, the system can easily 
survive for the age of the star. In fact, the plan- 
etary masses could be up to ~ 1.9 times bigger 
than their nominal values without violating stabil- 
ity constraints (Fig. fTTj) . This value is similar to the 
maximum planetary masses that can stably exist at 
the fi xed point of the 2:1 resonance (jBeauge et al.l 
l2003l . fig. 8). It will be very interesting to find a 
test of this hypothesis. The Laplace angle can also 
librate in this system, though this is not a require- 
ment for stability. 

Strong apsidal alignment, while unnecessary for system 
stability, allows planetary survival in a weaker 2:1 mean- 
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motion resonance. 

One final type of system architecture is, in princi- 
ple, consistent with stability of the HR 8799 system. 
The planets may be hierarchically spaced, with large ec- 
centricities and mutual inclinations, but inhabit phases 
of their orbits that look closely-packed at the moment. 
Such systems fit the current data but are finely tuned 
both in their orbital parameters and in the time at which 
we are viewing the system. 

It is possible that other stabilizing resonant configura- 
tions exist, but were missed because they occupy small 
regions of phase space. The 2:1 mean-motion resonance 
dominates our randomly-generated stable systems and 
would naturally yield the observations without fine tun- 
ing, so we consider it to be the most likely stabilizing 
mechanism. 

This study brings up several issues for future observa- 
tions of HR 8799, and directly-imaged multiplanet sys- 
tems in general, as follows. 

First, it serves as the first test of hot-start cooling mod- 
els for exoplanets. They barely pass the test if only 
planets d and c are in resonance, and they comfort- 
ably pass the test if all three planets are in resonance. 
We hope more detailed dynamical studies of this sys- 
tem will sharpen this test as more data are collected. If 
the doubly-resonant configuration can somehow be ver- 
ified, it would considerably weaken this test. We have 
not yet directly used dust observations as a dynamical 
constraint. The spectral energy distribution reveals a 
massive debris disk surrounding the planetary system, 
with an orbital radius of > 66 AU (.Williams & Andrews' 
l2006f ). Given that planet b is observed at Ob > 68 AU, 
we expect that the inner edge of the debris disk must 
be > 90 AU, and that future measurements and model- 
ing will find that orbital radius to be plausible and even 
preferred. Such a model could in turn serve as a comple- 
mentary test of planet b's mass, in analogy to the test of 
the mass of Fomalhaut b (jChiang et all 120091 ). Perhaps 
other directly-imaged systems will fortuitously arrange 
for complementary tests. 

Second, we found evidence of a mean motion resonance 
at very large orbital separations, much farther than those 
found by the ra dial velocity techn i que, of which ther e 
are many (e.g., iMarcv et all 1200 It iMavor et ani2004f ). 
Sometimes resonant identification is based on stabilit y 
arguments in those systems (e.g., ICorreia et al.l [20051 ) . 
as it is here. The most commonly invoked evolution- 
ary mechanism for trapping planets into resonance with 
one another is convergent migration in the protoplane- 
tary disk. The properties of migration in a disk with 
multiple massive planets deserve further investigation to 
determine the conditions under which convergent migra- 
tion and resonance capture are possible at the locations 
of the HR 8799 planets. We verified that if planets d and 
c were initially placed in circular orbits exterior to the 
2:1 resonance, they are stable to collisions or ejections 
for > 30 Myr, so getting into the resonance without first 



becoming unstable is not a problem in this case. One dif- 
ficulty with differential migration is that any additional 
migration, after the resonance is reached, efficiently in- 
creases the eccentricities. That this requires implausi- 
ble fine-tuning in the absence of eccentricity damping 
by the gas disk has been discussed for the 2:l-resonant 
GJ876 system (|Lee fc Pealell2002t ). In a trial integra- 
tion, we introduced a force to sim ulate outward migra- 
tion of the inner planet, following iLee fc Peald ([2002), 
with a timescale a/{da/dt) — 10^ yr and no eccentric- 
ity damping, starting from double-2:l resonance at the 
nominal masses (fig. [10], table |3|). The planetary ec- 
centricities rapidly increased and the system began scat- 
tering after a semi-major axis change of ^ 15%, illus- 
trating its fragility. We expect calculations of migration 
into resonance will be very interesting for this system. 
We also showed how perturbations by a third planet 
tend to disrupt a mean-motion resonance (when only 
the inner sub-system is in resonance). This mechanism 
adds to a growing list of ways to disrupt resonances 
among planets, including turbulent fluctuations in a pro- 
toplanetary disk ([A dams ct al. 2008), tidal dissipation 
(jTerquem fc Papaloizou ,2007,). a nd scat tering of plan- 
etesim als ( Murrav-Clav fc Chiang 2006; Morbidelli et al.l 
[200l . 

Third, it may seem surprising that dynamical stability 
arguments are needed to correctly solve the orbits of the 
first directly-imaged multiplanet system. However, this 
situation is also common for multipla net systems discov - 
ered by radial veloci ty. For instance, IVogt et aD ()2005l ) 
and iLee et al.l ([20061 ) have found that dynamical stabil- 
ity can constrain orbital parameters more tightly than 
radial-velocity data alone. Furthermore, planets that 
are discovered by direct imaging of their self-luminosity 
are biased to have high masses, making stability less as- 
sured. The bias of direct imaging towards large angular 
separation implies that very long orbital periods will be 
common for such discoveries, so we foresee many stabil- 
ity analyses predicated on only the sky-projected posi- 
tions and velocity vectors of planets. We hope this paper 
proves to be a useful example of how to conduct such an 
analysis. 
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Log-book of integrations: L Orbital Elements 
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0° 


180° 


0° 


0° 








±10-"* 


±10-4 


±0.01° 


±0.01° 


±0.01° 


±0.01° 




1.5 


0.0048 


67.91 


0.0 


I 


0° 




62.36" 


Fig. [IS 




0.0067 


37.97 


0.0 


0° 


0° 


0° 


316.99° 






0.0067 


24.44 


0.0 


i 


0° 




200.23° 



Note. — For each integration or suite of integrations listed, the three hnes specify the 
initial conditions for planets b, c, ari d d respectively. These are the values we input to the 
integrator Mercury ()Chaniberslll999[ ). 
(*) Mean anomaly. 

(**) In this series of integrations the planetary masses were scaled by various factors (x). 
To generate Figure ITl] small, random (Gaussian) components with the indicated standard 
deviation were added to the orbital elements for a statistical sample of the chaotic outcomes. 
(***) This series of integrations spanned a grid of three orientation parameters, i (common 
to planets b and d), fif,, and Qd- The values were i G [0°, 2°, 4°, 20°, 30°, 40°, 180°]; 
and Qd e [0°, 60°, 120°, 300°]. 
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TABLE 4 

Log-book of integrations: IL State Vectors 



lilt cgrcit ion 


A/* 


M„ 

ivip 


Xj\[ 










vz 




(M, 


si 




{ ATT^ 




(10 


-3 ATT ^-11^- 

AU day 


) 


^ nominal 


1.5 


0.0067 


60.16 


31.50 


0.0 


1.186 


-2.265 


0.0 


M A 




0.0095 


-25.90 


27.76 


0.0 


2.500 


2.332 


0.0 


Fig. [ni 




0.0095 


-8.45 


-22.93 


0.0 


-3.998 


1.473 


0.0 


ii3j 


1.5 


0.0067X 


60.16 


31.50 


0.0 


1.186 


-2.265 


0.0 


Fie: IgT' 




0.0095X 


—25.90 


27.76 


0.0 


2.496 


2.329 


0.0 






0.0095X 


-8.45 


-22.93 


0.0 


-3.996 


1.473 


0.0 




1.5 


0.00435 


60.16 


31.50 


0.0 


1.186 


-2.265 


0.0 


Fig. [7l8l 




0.0062 


-25.90 


27.76 


0.0 


2.496 


2.329 


0.0 






0.0062 


-8.45 


-22.93 


0.0 


-3.9967 


1.4737 


0.0 




1.60905 


0.0048 


60.18 


31.47 


27.21 


1.464 


-2.244 


0.4295 


Fig. [12] 




0.0067 


-25.94 


27.89 


-8.527 


2.398 


2.460 


1.178 






0.0067 


-8.601 


-22.90 


-5.311 


-3.577 


1.891 


-1.327 




1.54796 


0.0048 


57.54 


30.22 


9.765 


1.293 


-2.034 


0.2204 


Fig. [H 




0.0067 


-24.74 


26.61 


-23.15 


2.228 


2.327 


1.334 






0.0067 


-8.084 


-22.00 


-9.609 


-3.490 


1.518 


0.2634 



Note. — (*) In this series of integrations the planetary masses were scaled 
by various factors (x). Some reported integrations had only two planets: they 
were missing either the first planet listed (b) or the last planet listed (d). 
(**) In this series of integrations the innermost planet had a non-zero eccentric- 
ity. This was implemented by multiplying its initial velocity from the face-on, 
coplanar model by 7 = a/2 — 24.44AU/ad. Two series were performed, one with 
all three planets and one without the first planet listed (b). 



